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Abstract 

An explicit series solution is proposed for the inversion of the spherical mean 
Radon transform. Such an inversion is required in problems of thermo- and photo- 
acoustic tomography. Closed-form inversion formulae are currently known only for 
the case when the centers of the integration spheres lie on a sphere surrounding the 
support of the unknown function, or on certain unbounded surfaces. Our approach 
results in an explicit series solution for any closed measuring surface surrounding a 
| region for which the eigenfunctions of the Dirichlet Laplacian are explicitly known — 

such as, for example, cube, finite cylinder, half-sphere etc. In addition, we present a 
fast reconstruction algorithm applicable in the case when the detectors (the centers 
of the integration spheres) lie on a surface of a cube. This algorithm reconsrtucts 
3-D images thousands times faster than backprojection-type methods. 

> ' 

cN Introduction 

^ | The problem of image reconstruction in thermo-acoustic and photo-acoustic tomography 

is equivalent to recovering a function from a certain set of its spherical means [11, 13, 14, 
18, 19]. The process starts with object of interest being excited by a short electromagnetic 
pulse. This causes thermal expansion of the tissue, and generates an acoustic wave whose 
intensity is recorded by a set of detectors located outside the object. The intensity of the 
thermal expansion depends on the local properties of the tissue, and is of interest to a 
doctor, since abnormally high values of this function are indicative of a tumor. Under 
certain simplifying assumptions, the measurements can be related to the integrals of the 
expansion intensity over the spheres with the centers at the detectors' locations. The 
reconstruction of the local properties from these integrals is equivalent to the inversion of 
the spherical mean Radon transform. 

Some of the recent results on the injectivity of this transform as well as the correspond- 
ing range conditions can be found in [1-4,8,9]. In the present paper we concentrate on 
inversion formulae and algorithms for the solution of the reconstruction problem. Gener- 
ally, in such applications as photo- and thermo- acoustic tomography, the designer of the 
measuring system has a certain freedom of choice when selecting the detectors' positions 
(the centers of the integration spheres). Most of the known explicit solutions pertain to 
the spherical acquisition geometry, in other words to the configuration in which the de- 
tectors are located on a sphere surrounding the object. Such are the recently found series 
solutions [16-18] and backprojection-type formulae [8,10,12,19]. Explicit reconstruction 
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formulae are also known for such acquisition geometires as an infinite plane [7, 15, 19] and 
an infinite cylinder [19]. 

The spherical geometry is preferable to unbounded measuring surfaces since the latter 
have to be truncated in practice, which leads to errors in the reconstruction. However, 
there are compelling reasons to consider other non-spherical bounded measuring surfaces 
as well. For example, as shown in section 2, if the detectors are located on a surface of a 
cube surrounding the object of investigation, it is possible to design a fast algorithm that 
reconstructs the unknown 3-D function in a matter of seconds — as opposed to several 
hours required for the algorithms based on straightforward discretization of one of the 
3-D backprojection-type inversion formulae [8, 12, 19]. 

We thus present a series solution for the inversion of spherical mean Radon transform 
in the case when the centers of the integration spheres lie on a closed surface surrounding 
a bounded connected region in M n , n > 2. Our procedure requires knowledge of the 
eigenfunctions of the Dirichlet Laplacian defined on the region enclosed by the measuring 
surface. For many regions of practical interest such eigenfunctions are known explicitly. 
Among such regions in 3-D are, for example, a rectangle, a ball, an ellipsoid, a cylinder, 
a spherical shell. In addition, these eigenfunctions can be easily found for certain subsets 
of these bodies obtained by dissecting them along a plane of symmetry — for example for 
a half-ball, half-cylinder, certain triangular prisms and tetrahedra. Yet another example 
of regions with explicitly known eigenfunctions is given by the crystallographic domains 
(see [5,6] for details). A generalization of this approach to a general connected region 
is possible if one computes the eigenfunctions of the Dirichlet Laplacian numerically. In 
this case, however, the reconstruction algorithm is likely to be rather expensive from the 
computational point of view. 

The proof of the range theorem in [1] involves implicitly a reconstruction procedure 
also based on eigenfunction expansions. Unlike the present method, that procedure would 
involve division of analytic functions that have countable number of zeros. While the range 
theorem guarantees cancellation of these zeros when the data are in the range of the direct 
transform, a stable numerical implementation of such division would be complicated if 
not impossible. (A similar problem arises with the series solution of [16] that involves 
division of certain computed quantities by Bessel functions.) The technique we present 
below does not require such divisions. 

Section 1 contains a general description of the present method. The efficiency of 
numerical realization of this technique depends, in particular, on the availability of fast 
algorithms for the summation of the arising eigenfunction expansions. In the simplest 
case of a cubic (or rectangular) measuring surface such an algorithm is the 3-D Fast Sine 
Fourier transform. This allows us to design a very efficient reconstruction algorithm for 
this particular measuring configuration, as discussed in section 2. Finally, in section 3 
we investigate an interesting property that seems to be exclusive to the series solutions 
presented in this paper. Namely, this technique will produce a theoretically exact image 
within the region enclosed by the measuring surface even if there are sources outside that 
region. This property can prove to be useful for reducing the sensitivity of the measuring 
system to external noise. Such a noise cancellation will occur, however, only if all the 
measurements are performed simultaneously by a fixed set of detectors; a synthesized 
measuring surface will not exhibit this phenomenon. 
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1 Series solution 



Suppose that Cq function /(x), x eW 1 , n > 2 is compactly supported within the 
bounded connected open region f2 with boundary <9f2. Our goal is to reconstruct /(x) 
from its projections g(z,r) defined as the integrals of /(x) over the spheres of radius r 
centered at z: 



<?(z,r) = y /(z + rfV" 1 ^*), 



jn-l 



where § n is the unit sphere in R n , t is a unit vector, and ds is the normalized measure 
in MJ 1 . Projections are assumed to be known for all z e <9f2, < r < diam(fi) (integrals 
for r > diam(Q) automatically equal zero, since the corresponding integration spheres do 
not intersect the support of the function). 

Suppose A^, M m (x) are the eigenvalues and normalized eigenfunctions of the Dirichlet 
Laplacian —A on Q with zero boundary conditions, i.e. 

Aw m (x) + \ 2 rn u m {x) = o, xeft, ficr, (1) 
■u m (x) = 0, x G 

||«m||2 = / \u m M\ 2 d*L = 1. 



We would like to reconstruct function /(x) from the known values of its spherical integrals 
g(z,r) with the centers on dQ: 



g(z,r)= J f(z + rs)r n 1 ds, zGffi. 



Urn. \ X ) 



We notice that u m (x) is the solution of the Dirichlet problem for the Helmholtz equation 
with zero boundary conditions and the wave number A m , and thus it admits the Helmholtz 
representation 

/ $ Am (]x-z|)— u m {z)ds{z) xGfl, (2) 
Jan on 

where <3>A m (l x — z|) is a free-space rotationally invariant Green's function of the Helmholtz 
equation (1). 

Our approach is based on the fact that eigenfunctions {■u m (x)}^° form an orthonormal 
basis in L 2 (f2). Therefore /(x) can be represented by the series 

oo 

/( x ) = a « M ™( x ) ( 3 ) 

m=0 

with 

a 



= / ii m (x)/(x)dx. (4) 
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Since /(x) is Qj, series (3) converges pointwise. The reconstruction formula will result if 
we substitute representation (2) into (4) and interchange the order of integrations 



= / u m (x)/(x)dx 
Jn 

= In {Ian $Am(|X ~ z ^^ Um ^ ds ^) f ^ d * 
= J Q[ $ A ra (|x - z|)/(x)dx) ^-u m (z)ds(z) (5) 
d 

I{z,X m )—u m (z)ds(z), (6) 
an ° n 

/(z,A m ) = / $ Am (|x-z|)/(x)dx. 

The change of the integration order is justified by the continuity of eigenf unctions u m (x). 
Function /(z, \ m ) is easily computed from the projections 

J(z,A m ) = f $ Am (|x-z|)/(x)rfx= f g(z,r)<S> x Jr)dr, 
Jn J 



where 



R+ 



and with Fourier coefficients a m now known, /(x) is reconstructed by summing series (3). 
If desired, this solution can be re-written in the form of a backprojection-type formula: 

/( X ) = Y] «m M m( X ) = / «m$A m (|x - z\)—U m (z) ds(z) 

m =o Jdn \ m =o dn J 

= [ h(z,\x-z\)ds(z), (7) 
Jan 



where 



°° d 
h{z,t) = ^a m $ Am (t)— w m (z), 



dn 

m=0 



and coefficients a m are computed using equation (6). In the above formula equation (7) 
is clearly a backprojection operator, and (8) is a filtration. However, the latter operator 
is now represented by a series rather than by a closed form expression. Moreover, this 
operator is not local in z, unlike the filtration operator of the known closed-form explicit 
inversion formulae [8, 10, 12, 19]. 

Finally, we notice that if function /(x) is not smooth but rather belongs to L 2 (Q), our 
reconstruction formulae are still valid if equation (3) is understood in the L 2 sense. 



2 A fast algorithm for the cubic measurement surface 
in 3D 

A cube is the simplest of the regions whose eigenfunctions of the Dirichlet Laplacian are 
known explicitly; they are products of sine functions. In the present section we exploit 



4 



the simple structure of these eigenfunctions to develop a fast reconstruction algorithm 
applicable in the case when the detectors are located on a surface of a cube (a general- 
ization to a rectangular case is straightforward). Such a measuring surface can be either 
sampled by regular detectors or synthesized from measurements made by interferometric 
line detectors as discussed in the Introduction. 

Let the sought function /(x) be supported within the cube Q = [0, R] x [0, R] x [0, R}. 
We will index the normalized eigenfunctions w m (x) and eigenvalues A m of the Dirichlet 
Laplacian on this region using vector m = (mi,m 2 ,m 3 ), mi,m 2 ,m3 G N : 



Uy 



Am 



o imiiXi nm 2 x 2 nm 3 x 3 
—r sin — - — sin — - — sin — - — . 
R 3 R R R 

2i |2 
7T m . 



Cube Q has six faces Sfli, i — 1, 6 : 

5Qi = {x|xi = R, < x 2 < R, < x 3 < R, }, 

5il 2 = {x.\xx = 0,0 < x 2 < R,0 < x 3 < R,}, 

5Q 3 = {x\x 2 = R, < xi < R, < x 3 < R, }, 

5fi 4 = {x\x 2 = 0,0 < x 1 < R, < x 3 < R, }, 

5Q 5 = {x|z 3 = R, < xi < R, < x 2 < R, }, 

SQq = {x|x 3 = 0,0 < x 1 < R, < x 2 < R, }. 

The values of the normal derivatives ^ii m (x) of the eigenfunctions on the boundary are 
equal to certain products of sine functions: 



d_ 
dn 



(-1) 



87rmi • Trm 2 x 2 ■ Tvm 3 x 3 

R i sin — R — sin — R — 

iii i 8tttoi g j n Trm 2 x 2 g j n Tvm 3 x 3 
8ttto 2 ■ wrnixi ■ wm 3 x 3 

Ri sin R sm R 

v L) Ri Sill R sm H 

87tto 3 • wm 1 x 1 ■ Trm 2 x 2 

R A R R 

8nm 3 • nmixi ■ ■wm 2 x 2 

IF™ R bU1 R 



-1 



\m 3 



x e5Qi 
x e5tt 2 

x G#fi 3 

X G(%7 4 
X G^5 
X G5^6 



(9) 



As in section 1, in order to reconstruct /(x) we recover Fourier coefficients a x 

f d 

«m = / I(z,\ m )—U m (z)ds(z) 

Jan ° n 

6 f d 
= / /(z, A m )— w m (z)rfs(z), 



where 



7(z,A) 



#(z, r)&\(r)dr. 



(10) 



(11) 



If we choose the Green's function 3>\(t) in the form 

cos At 



Ant 



equation (11) can be re-written in the form of the Cosine Fourier transform as follows 





~g(z,r)~ 


4tt / 


r 



V3R_ 

nil. Til 

cos Xrdr. (12) 

As before, when coefficients a m have been found function /(x) is obtained by summing 
the Fourier series 

/(x) = ^ a m«m(x). (13) 
meN 3 

The above formulae are just a particular case of the inversion technique presented 
in the previous section. They yield theoretically exact reconstruction if the effects of 
discretization are neglected. However, in the practical computation only limited range 
of frequencies < A < \ N yi mst C an be recovered from finitely sampled (in r) projections 
g(z,r) using equation (12). Therefore series (13) has to be truncated. The Gibbs phe- 
nomenon resulting from such a truncation can be reduced by application of a filter r/(A m ), 
so that instead of the previous equation the following formula will be used to reconstruct 
an approximation to /(x) : 

/(x) w ^ « m ^(A m )M m (x). (14) 

meN s ,\\ m \<\ Nyqulst 

The whole reconstruction procedure can be accelerated by utilizing the Fast Cosine 
Fourier transform to compute (12), the 3-D Fast Sine Fourier transform to sum series 
(13), and the 2-D Fast Sine transform to evaluate the six integrals in equation (10). 
However, there is one obstacle for implementing this plan. The integrals in (10) need to 
be computed for different values of A m , and there are too many of these values to make 
the algorithm efficient. The work-around for this problem is to evaluate these integrals for 
a set of uniformly distributed values A; = I A A, / = 0, 1, 2, ... and then to find the needed 
values for each of A m by interpolation. Such an interpolation in the spectral parameter A 
requires careful selection of discretization steps and interpolation techniques. The details 
of our implementation are presented below. 

Suppose function /(x) is to be reconstructed on n x n x n Cartesian grid, and the 
detectors are located at the nodes of 2-D n x n Cartesian grids defined on the faces of 
the cube (values at the edges of the cube will not be needed). We will assume that the 
discretization step of measurements (in r) is approximately the same as the step of the 
Cartesian grids. Then the number of samples n\ in one projection g(z,r) approximately 
equals \^3n. Depending on a type of the Fast Cosine Fourier transform algorithm used 
to compute (12), the projections will have to be padded with zeros to make the total 
number of samples n 2 either a power of 2, or a product of small prime numbers. Thus, 
rii < n 2 < 2n\. The Nyquist frequency \ N yi mst corresponding to such discretization equals 
n(ni — l)/D, where diameter D = diam(Q) = y/3R. With these parameters in mind we 
summarize the five steps of the fast reconstruction algorithm. 

Step 1. The first step is to compute a discrete version of (11) using the Fast (discrete) 
Cosine Fourier transform of length n 2 . This will produce values /(z, A;) for the frequencies 
A; = l\ Nyqmst / {n 2 — 1), I = 0, 1, ..,n 2 — 1. Importantly, as long as n\ < n 2 , the step of 
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discretization of I(z, \) in A is small enough to make it possible to approximately recover 
values of I(z, A) (or a linear function of I(z, A)) for A 7^ A/ by interpolation. 
Step 2. For each value of A/ compute integrals in the form 




(15) 



for integer all integer nii, m k by means of the 2-D Fast Sine Fourier transform. 
Step 3. Compute approximate values of integrals 




by interpolating values obtained on step 2 (equation 15). Some care should be taken 
to guarantee accuracy of computations on this step. It is well known that a low or- 
der interpolation in spectral parameter can lead to a suboptimal reconstruction, due to 
highly oscillatory nature of the Fourier transformant. An example and analysis of this 
phenomenon can be found in [15]. However, the Fourier transform of a finitely supported 
function is an analytic function of the spectral parameter, even if the function itself is 
known imprecisely. Therefore, higher order polynomial interpolation is applicable and 
does produce good results in this case. In our numerical experiments we observed that, 
indeed, the linear interpolation in A yields rather inaccurate reconstruction. The increase 
in the order of the polynomial interpolation significantly improves the image; if the 6- 
th order Lagrange interpolation is utilized, the interpolation error is dominated by the 
discretization errors and further increase in the accuracy of interpolation is not needed. 

Step 4. Use values computed on step 3 to calculate A m by combining equations (9) 
and (10). 

Step 5. Using the 3-D Fast Fourier Sine transform to implement (14), compute values 
of /(x) at the nodes of the 3-D Cartesian grid. 

A simple computation shows that the number of floating point operations implemented 
on each steps of the algorithm is 0(n z \ogn) for steps 1,2, and 5, and 0(n 3 ) on steps 3 
and 4, resulting in a total 0(n 3 \ogn) operation count for this technique. This has to 
be compared with 0(n 5 ) operation count required by a backprojection step of a method 
resulting from a straightforward discretization of any of the explicit inversion formulae[8, 
12,19]. (The latter estimate assumes that the reconstruction is done on n x n x n Cartesian 
grid from n 2 detector positions.) 

In what follows we present a numerical example illustrating the work of our algorithm. 
A function is reconstructed within the cube [0, 1] x [0, 1] x [0, 1]. The dimension of the 
grids were defined by the values of parameters n = 129, n\ = 223, and n<i = 256. As a 
filter we used the cosine window function 



We utilized the same phantom as in [12] to facilitate the comparison of the present results 
with the images reconstructed in the former work by application of discretized explicit 
inversion formulae. The phantom consists of eight characteristic functions of the balls 




A < \ N vi uist 
A > \ N yi uist 
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with radii ranging from 0.06 to 0.13, whose centers lie in the plane £3 = 0. The cross 
section of the phantom by the latter plane is shown in Figure 1(a). Figure 1(b) shows the 
central cross section (x 3 = 0) of the reconstruction from the exact data. 




(a) (b) (c) 

Figure 1: Numerical example: (a) the phantom, (b) reconstruction from the exact data, 
and (c) reconstruction from noisy data 

In order to evaluate the sensitivity of the algorithm to noise in data, the imprecise 
measurements were modeled by adding to the projections normally distributed noise with 
the intensity 15% of the signal (in L2-norm). In this experiment the values of the recon- 
structed function were set to zero outside of the cube [0.05, 0.95] x [0.05, 0.95] x [0.05, 0.95], 
since the reconstruction from noisy data is unstable at the locations close to the detectors, 
due to the singular nature of the Green's function. Such sensitivity is natural, and is an 
issue for other reconstruction techniques as well. For example, the slices of 3-D images 
obtained in [12] were computed within the unit ball from the detectors located on a sphere 
of radius 1.1; the reconstruction in the close vicinity of the detectors was also avoided. 
The gray scale representation of the reconstructed function is shown in Figure 1(c); in 
Figure 2 we demonstrate the surface plot of the same function. 

We believe that the quality of the reconstructed images in all of the above experiments 
is as good as of those resulting from explicit inversion formulae (see [12]). Similarly to 
those formulae the present technique demonstrates high stability of the algorithm to the 
perturbations of the data. On the other hand, the computation time for the present 
algorithm in the above experiment varied from 7 to 8 seconds on an AMD workstation 
with a 2GHz processor; function values at about 2 million points were reconstructed from 
about 97 thousand projections. For comparison, the reconstruction reported in [12] from 
33000 projections at about a million grid points by the fastest of our implementations of 
the explicit inversion formula took about 48 minutes. If the latter algorithm is used to 
process the same number of the projections and grid points as in the present example (97 
thousands and 2 million respectively), the computation time increases to about 7 hours. 
It is fair to say that the fast algorithm is thousands time faster than the straightforward 
discretization of any of the backprojection-type formulae. 
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Figure 2: Reconstruction from noisy data; surface plot 

3 Reconstruction in the presence of exterior sources 

The series solution described above has an interesting property not possessed (to the best 
of our knowledge) by any other currently known explicit reconstruction technique. Let us 
consider a slightly more general problem. Suppose that region Q is a proper subset of a 
larger region Qi (Q C and that a L 2 function F is defined on f^. We will denote the 
restriction of F on Q by /, i.e. 

" F(x) , x6ll 
, x G R n \tt ' 



/(x) = 



We would like to reconstruct /(x) from the integrals g(z,r) of F over spheres with the 
centers on dVt: 



g(z,r)= J F(z + rs)r n_1 ds, zGffi. 



Unlike in the previously considered problem, now the centers of the integration spheres 
are lying on a surface contained within the support Qi of the function F. While we are 
still trying to reconstruct the restriction / of F to Q, the integrals we know are those of 
F and not of /. 

It turns out that the solution to this problem is still given by formulae (3) and (6) 
(or equivalently by (7), (8), and (6)). Indeed, if we extend functions u m (x) by to M. n , 
formula (2) holds for all x G K n , and (3) remains unchanged. In formula (5) / can be 
replaced by F as follows: 



u m (x)/(x)dx = / u m (x)F(x)dx 
Jn 

= J (j^Aj|x-z|)F(x)dx) |- Mm (z)rf S (z), 

and the inner integral can be computed from projections as before: 

/ $A m (|x-z|)F(x)dx= / g(z,r)® Xm (r)dr. 
Jn J 
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(a) 



(b) 



Figure 3: Reconstruction in the presence of exterior sources (a) phantom; the white line 
shows location of the detectors (b) reconstructed image. 

By combining the two above equations we again arrive at the formula (6). 

This interesting property can be illustrated by a numerical example. We consider the 
same phantom as in the previous section. This time, however, the detectors are located 
on the surface of the cube Q = [0.235,0.765] x [0.235,0.765] x [0.235,0.765]. Location 
of the detectors is shown in Figure 3(a) by a white line. The integrals were computed 
over full spheres, and the reconstruction was conducted, as before on the grid of size 
129 x 129 x 129 within Q. The result is presented in Figure 3(b), and, as a surface plot, 
in Figure 4. On the latter figure one can notice shallow troughs resulting from imperfect 
resolution of sharp edges of exterior sources by a finite number of detectors. The depths 
of these troughs, however, does not exceed 6% of the maximum of the original function. 




Figure 4: Reconstruction in the presence of exterior sources; surface plot. 



10 



4 Acknowledgments 



The author would like to thank P. Kuchment for fruitful discussions and numerous helpful 
comments. 

This work was partially supported by the NSF/DMS grant NSF-0312292 and by the 
DOE grant DE-FG02-03ER25577. 

References 

[1] M. Agranovsky, P. Kuchment, and E. T. Quinto, Range descriptions for the spherical 
mean Radon transform, preprint 2006, arXiv: math. AP/0606314. 

[2] M. L. Agranovsky and E. T. Quinto, Injectivity sets for the Radon transform over 
circles and complete systems of radial functions, Journal Of Functional Analysis, 
139, no. 2 (1996) pp. 383-414. 

[3] G. Ambartsoumian and P. Kuchment, On the injectivity of the circular Radon trans- 
form arising in thermoacoustic tomography, Inverse Problems 21 (2005), pp. 473-485. 

[4] G. Ambartsoumian and P. Kuchment, A range description for the planar circular 
Radon transform, SIAM J. Math. Anal. 38, no. 2 (2006) pp. 681-692. 

[5] P. Berard, Spectres et Groupes Cristallographiques, C. R. Acad. Sci. Paris A-B 
288, no. 23 (1979) pp. A1059-A1060. 

[6] P. Berard and G. Besson, Spectres et Groupes Cristallographiques II: Domaines 
Spheriques, Ann. Institut Fourier, 30, no. 3 (1980) pp. 237-248. 

[7] J. A. Fawcett, Inversion of iV-dimensional spherical averages, SIAM J. Appl. Math. 
45, no. 2, (1985) pp. 336-341. 

[8] D. Finch, Rakesh, and S. Patch, Determining a function from its mean values over a 
family of spheres, SIAM J. Math. Anal. 35 no. 5 (2004), 1213-1240. 

[9] D. Finch and Rakesh, The range of the spherical mean value operator for functions 
supported in a ball, Inverse Problems 22 (2006), pp. 923-938. 

[10] D. Finch, M. Haltmeier, and Rakesh, Inversion of spherical means and the wave 
equation in even dimensions, preprint 2006. 

[11] P. Kuchment, Generalized Transforms of Radon Type and Their Applications, in G. 
Olafsson and E. T. Quinto (Editors), The Radon Transform, Inverse Problems, and 
Tomography, Proc. Symp. Appl. Math. v. 63, AMS, Providence, RI 2006, pp.67 - 91. 

[12] L. A. Kunyansky, Explicit inversion formulae for the spherical mean Radon transform, 
to appear in Inverse Problems. 

[13] R. A. Kruger, P. Liu, Y. R. Fang, and C. R. Appledorn, Photoacoustic ultrasound 
(PAUS) reconstruction tomography, Med. Phys. 22, (1995), pp. 1605-1609. 



11 



[14] R. A. Kruger, D. R. Reinecke, and G. A.Kruger, GA Thermoacoustic computed 
tomography - technical considerations, Med. Phys. 26, no. 9, (1999) pp. 1832-1837. 

[15] F. Natterer and F. Wubbeling, Mathematical Methods in Image Reconstruction, 
Monographs on Mathematical Modeling and Computation, 5, SIAM, Philadelphia, 
PA 2001. 

[16] S. J. Norton, Reconstruction of a two-dimensional reflecting medium over a circular 
domain: Exact solution, J. Acoust. Soc. Amer. 67 (1980), pp. 1266-1273. 

[17] S. J. Norton and M. Linzer, Ultrasonic Reflectivity Imaging in Three Dimensions - 
Exact Inverse Scattering Solutions for Plane, Cylindrical, and Spherical Apertures, 
IEEE Trans. Biomed. Eng. 28, no. 2 (1981) pp. 202-220. 

[18] M. Xu and L.-H. V. Wang, Time-domain reconstruction for thermoacoustic tomog- 
raphy in a spherical geometry, IEEE Trans. Med. Imag. 21, (2002), pp. 814-822. 

[19] M. Xu, L. V. Wang, Universal back- projection algorithm for photoacoustic computed 
tomography, Phys Review E 71 (2005), p. 016706. 



12 



